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1 Introduction 

This paper presents and derives the interrelations between survival analysis and master 
equation. Both have important applications in the social sciences and other scientific fields 
treating stochastic systems. However, since they focus on different aspects of modeling, it 
is not yet generally known that they are closely related to each other. 

Survival analysis deals with modeling the transitions between succeeding states of a 
system. Questions related with this are the timing, spacing, and sequencing of the states 
of a time series. Survival analysis tries to fit and understand the distribution of these 
quantities in terms of the functional form of the hazard rates which are responsible for 
the investigated transitions. The parameters specifying the concrete functional form of 
the hazard rates are normally estimated from empirical data by means of the maximum- 
ox partial-likelihood method. 

Once the hazard rates are empirically known or hypothetically specified, one can carry 
out microsimulations of corresponding time series by means of the Monte-Carlo simulation 
method. This allows the prognosis and investigation of the characteristics and frequency 
of time series. 

However, if one is interested in cross-sectional data connected with the stochastic 
process under consideration, one needs to know the temporal evolution of the distribution 
of states. This can be obtained by simulation of the associated master equation, which 
only presupposes that the initial distribution and the hazard rates are given. In other 
words: The master equation is suitable for the calculation and prognosis of cross-sectional 
data which are related with the longitudinal life-table data used for survival analysis. 

In addition, some new formulas are introduced which allow the determination of path- 
related (i.e. longitudinal) quantities like the occurence probability, the occurence time 
distribution, or the effective cumulative life-time distribution of a certain sequencing of 
states (path). These can be efficiently evaluated with a recently developed simulation tool 
(EPIS) which also provides a new solution method of the master equation (the contracted 
path integral solution). In contrast, a calculation on the basis of time-series data would 
require an extreme computational effort. 

The effective cumulative life-time distribution facilitates the formulation of a hid- 
den state concept of behavioral changes which allows an interpretation of the respective 
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time-dependence of hazard rates. Hidden states represent states which are either not phe- 
nomenological distinguishable from other states, not externally measurable, or simply not 
detected. They could, for example, reflect the psychological and mental stages preceding 
a concrete action, the individual predisposition, or the internal attitude towards a certain 
behavior. 

2 Survival Analysis 
2.1 Preliminaries 

Survival analysis deals with the investigation of the timing, spacing, and sequenc- 
ing of longitudinal time series (panel studies), e.g. life-table data (Blossfeld et al. 
1986; Kalbfleisch/Prentice 1980; Elandt- Johnson/ Johnson 1980; Diekmann/Mitter 1984; 
Tuma/Hannan 1984; Cox/Oakes 1984; Lancaster 1990; Courgeau/Lelievre 1992). These 
time series have the general form 

Y(t) = Y n for T n <t<T n+u (1) 

where T n is the variable of the waiting time after which the state Y n is occupied. The state 
Y n remains occupied for a time period T n < t < T n+1 ending directly before the next event 
(F n+ i,T n+1 ). Two succeeding events (Y n __i,T n _i) and (Y n ,T n ) determine the nth episode 
((F n _i,T n _ 1 ), (Y n ,T n )) of a time series. 

The state Y n stems out of a set Ai n of possible outcomes which are often called risks. 
Without loss of generality^] we can assume that this set does not depend on the respective 
episode n which implies M. n = Ai. Furthermore, M. is normally a discrete set so that 

Y n eM = {l,2,...,i,...,N}. (2) 

The elements i of this set represent the possible outcomes. Examples for Ai are: 

M. = {single, married, widowed, divorced} (3) 

or 

At = {school, army, training, studies, employment, unemployment, retirement} . (4) 

Of course, the sequence of waiting times T n = t n ('timing ') and often also the sequence 
of states Y n = i n ('sequencing') vary individually so that the actual time series 

Y a (t) = i n for t n <t < t n+1 (5) 

can only be described by a stochastic process 

(Y,T):={(Y n ,T n ):n = l,2,...}. (6) 

The subscript a G {1,2, ... ,A} distinguishes the different individuals, systems, or real- 
izations to which the respective time series belong. 

1 One can simply define M. := {JM n , where "(J" symbolizes the union of sets. If, for a given event 
history H n -\, no transitions take place to state j € M. one must only set the corresponding transition 
rate equal to zero, e.g. for j = i n -\. 
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Now, the spacing of a time series can be denned as the sequence of life times (survival 
times, failure times) 

V n :=T n -T n _ x , fc = l,2,.... (7) 

Additionally, we introduce the so-called event history of a state Y n = i n which has been 
occupied at time T n = t n by 

H n -i := {to, z'o; t%, i%, Xi; . . . ; t n _i, i n -i, x n _i} . (8) 

Here, the vector of covariates = (xk,ii ■ ■ ■ i x k,M) comprises the different factors Xk, m 
like education or sex which may influence the transition (i.e. change) from the (k — l)st 
state Yfe_i = Zfc_i to the kth state Y k = i k . 

The following figure illustrates the terms which were introduced above: 
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Figure 1: Example illustrating a time series with N = 4 competing risks and 3 episodes which 
are given by the sequence of events (1, to) — * (3,ti) — > (2,^2) — ► (4,^). 



2.2 Central Concepts 

We now come to the definition of some central concepts of survival analysis. In view of 
the following discussion, we restrict ourselves to the special case that 

1. only the last state Y n -\ = i n -i of the event history if n _i has an influence on the 
nth transition (Markov case), 

2. the transitions and the vector of covariates independent of the respective 
episode n. 

The second assumption implies that the vector of covariates x is time-independent during 
the survey. This means that each individual or system a is characterized by a fixed value 
of x and that different vectors of covariates distinguish different cohorts (subpopulations). 

Then, the so-called hazard rates (transition rates) X x (t;j\i) of subpopulation x are 
defined as the probability P x (T n < t + At, Y n = j\T n > t, Y n -i = z)Q per unit time At > 

2 Instead of P x (T n < t + At,Y n = j\T n > t,Y n _i = i) one sometimes writes P x (t < T n < t + 
At, Y n = j\T n > t,Y n -i = i). However, the expression t < T n < t + At is unnecessary complicated since 
t < T n is already presupposed by the condition T n > t in the second part of the argument of P x . 
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to change into state Y n = j up to time t + At on the conditions that the nth transition 
did not happen before time t and the preceding state was Y n -\ = i' 

A x (t; := ton ^^x(T n <t + At,Y n = j\T n > t, Y n ^ = i) . (9) 

Moreover, the survivor function (survival function) S™(t\i) of cohort x is defined by 
the probability P x (T n > £|l^_i = z) that the nth transition does not take place before 
time T n , given that the preceding state was Y n _i = i: 

S%(t\i) := P x (T n > t|y n _! = z) . (10) 

One can find the relation 

SJ(t|i) = 1 - F-(t\i) (11) 

with regard to the cumulative distribution function (life-time distribution, failure distri- 
bution, duration distribution function) 

F?(t\i) := P x (T n < t|Y n _! = z) (12) 

of subpopulation x. The latter describes the probability that the nth transition happens 
before time t, given that the preceding state was Y n _i = i. If we are interested in the 
probability that the nth transition not only happens before time t but additionally leads 
to state Y„ — j ^ i, we need the quantity 

F2(t;j\i) := P x (T n < t,Y n =j\Y n _ 1 = i) (13) 

which fulfils 

N 

Y,K&M = p x(T„ < t|Y n -i = i) = F?(t\i) . (14) 

3=1 

Finally we define the probability density function (failure time (sub) density function 
for failure type j) by 

f£{t;j\i) := Km Q -^P x (t<T n <t + At,Y n = j\Y^ l = i) 

P x (T n <t + At,Y n = j\Y n -! = z) - P x (T n <t,Y n = i|y n _! = i) 

= hm 

At^O At 

= J t K ( t;m- (15) 
(|) to ([l]]), an d (p~5|) imply the relations 

T,mm = j/mi) = -pzm (i6) 



3=1 



and 



\ x {t;j\i)££{t\i) = X x (t;j\i)P x (T n > t\Y n _ x = i) 

= lim -J-P x (T n < t + At, Y n = j\T n > t, y_! = z)P x (T n > t|Y n _ a = i) 

At ^0 /\t, 

= hm ^Pc(Pn < t + At, Y n = j, T n > £|y n _! = z) 

= a?) 
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Equations fll6|) and (fPTD yield the differential equation 
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which has the important solution 
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dt' A x (t';j|z) 
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with the overall departure rate ( overall hazard rate, overall failure rate ) 



(19) 
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(20) 



That is, survivor functions are determined by exponential relations which only depend on 
the hazard rates. 



3 The Master Equation Technique 

3.1 Derivation of the Master Equation 

We will now derive a system of differential equations for the evolution of the probability 
distribution P x (j,t) of states j with time t which is associated with the above described 
stochastic process. (Each subpopulation x obeys its own system of equations.) For this 
purpose we apply two relations from probability theory: 
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(21) 
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and 
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P x 0,t') = "£P x {j,t'\i,t)P x {i,t) . 



:22i 
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Here, P x (j, t) denotes the probability that we have state j at time t! given that we 



had state % at time t. With (glj) and fl22[) we get 
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Due to 



hm ^Px(j, t + At|z, t) = lhn ^Px(Tn < * + At, Y n = j\T n > t, Y n _ x = i) (24) 



and (0) we finally obtain the differential equation 



d N 

-P x (j,t) = £ [A x (t; j\i)P x (i,t) - X x (t;i\j)P x (j,t)] (25) 



which is called the master equation (Haken 1983; Weidlich/Haag 1983; Helbing 1995). Chi- 
ang (1968: 116ff), Lancaster (1990: pp. 109ff), and Courgeau/Lelievre (1992: 40) presented 
related considerations for continuous-time Markov processes leading to the Kolmogorov 
differential equation 

d N N 

Px{j,t\i ,t ) = ^T / X x (t;j\i)P x (i,t\i ,to) with A x (t; j\j) := - ^ A x (i; i\j) . (26) 



dt 



i=l «=i 



As known from the theory of stochastic processes, the master equation can be obtained 
from ( |2"B] ) by multiplication with the initial distribution of states P x (io, to) and subsequent 
summation over the initial states io because of P x (j,t\io,t )P x (i ,t ) = P x (j,t;i ,t ) = 
P x (Mo|j,t)P x (j,t)and©. 

According to the master equation, the temporal change of the probability P x (j,t) to 
have state j at time t is given by the sum of the effective transition rates X x (t;j <— i) 
from other states i to state j minus the sum of the effective transition rates X x (t; i <— j) 
from state j to other states i. The effective transition rate A x (t; j <— i) := X x (t; j\i)P x (i,t) 
from state i to state j is obviously the product of the transition rate A x (t; j\i) given that 
the indiviual or system under consideration is in state i times the probability P x (i, t) that 
he/she/it is actually in state i. 

In the following we will assume the special case that we have the dependence 

X x {t;j\i) := X° x {t)w x {j\i) (27) 

of the transition rates. Often one even restricts oneself to the proportional hazard rate 
model (proportional hazards model) 

X x {t;j\i) := A (t)exp(/3 iiX ) (28) 

so that A x (t) = Ao(t) with the baseline hazard function Xo(t) and w x (j\i) = exp(/3 ? jx). The 
optimal parameter vectors (3j i can be estimated from life-table data via the partial likeli- 
hood method proposed by Cox (1975) (see also Blossfeld et al. 1986; Kalbfleisch/Prentice 
1980; Elandt- Johnson/ Johnson 1980; Diekmann/Mitter 1984; Tuma/Hannan 1984; 
Cox/Oakes 1984; Lancaster 1990; Courgeau/Lelievre 1992). 

For reasons of simplicity we utilize relation (|27l) to introduce the subpopulation-specific 
times 

C dT 
T x (t) := / dt'X° x (t') corresponding to — - = A°(t) . (29) 

J (JjL 

to 

This implies, for example, 

pMt) - ±PJj,T M ] = (±PMrS) ^ - (skCi.rj) m . (30) 
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As a consequence, we can bring the master equation into the form 

d N 

-i— P*(j,T x ) = hxOIO-Px^Tx) -w x (i\j)P x (j,r x )] (31) 

" r x i=l 

with time-independent transition rates 

A x (r x ; = w x (j\i) . (32) 
Analogous simplifications are found for the other quantities: 

TV 

K(r x \i) = w x (i) := w x(jK) , (33) 

3 = 1 

N 

S:(r x \i) = exp [-w x (i)(r - rx>B _i)] = 1 - i£(r x |i) = 1 - £ F x "(r x ; , (34) 

and 



3 = 1 

(3*i) 



K(T x ;j\i)=w x (j\i)S:(T x \i) = ^F^(r x ;j\i), (35) 

dr x 

where we have used the convention r xn := r x (t n ). Since the subscript x for the subpopu- 
lation (cohort) is arbitrary but fixed (time-independent) we will omit it in the following 
which makes the mathematical formulas easier to read. 



3.2 Simulations with the Master Equation 

The master equation has proved to be a very powerful and quite flexible tool for the 
description and simulation of stochastically behaving systems which are subject to in- 
herent or external random influences ('fluctuations'). It has got numerous applications 
in physics, chemistry, biology, economics, and the social sciences (Weidlich/Haag 1983; 
Troitzsch 1990; Weidlich 1991; Helbing 1995). The master equation can normally not be 
analytically solved. However, since the master equation has the form of a linear system of 
ordinary differential equations, it can be easily simulated by means of the usual numerical 
integration algorithms (Press et al. 1992: Chap. 16). 

With today's computers these simulations are very fast, so that even systems with some 
dozens of variables can be solved within a few minutes. Therefore, the master equation 
technique facilitates an efficient evaluation of all cross-sectional quantities related with 
stochastic processes. This includes the temporal course of the distribution of states P(j, r), 
its maxima, means value, variance, etc. Whereas the mean value represents the average 
behavior of a huge number of time series, the course of the most probable time series 
will often (but not always) be close to the maxima of the distribution of states, and the 
variance is a measure for the grade of variation between different time series. In the case 
of a multimodal distribution of states, the mean value can considerably differ from the 
most probable time series. 

In summary, simulations of the master equation are suitable for scenario techniques, 
since they allow 
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1. to investigate and compare the properties of different conceivable stochastic models 
(which are normally related to different functional forms of the hazard rates), 

2. to evaluate the implications of parameter changes (which may correspond to con- 
sidered or planned modifications of legal regulations or social conditions), 

3. to make prognoses of the range of probable future behaviors of a stochastic system. 



4 Simulation of Path- Related Quantities 

4.1 Microsimulations with the Monte-Carlo Technique 

Although the master equation technique is very versatile, it does not allow investigations 
which are directly related with the longitudinal time series. However, the simulation of 
these is desirable for a number of reasons: 

1. The single time series give an impression of the possible time-dependent behavior 
of individual systems. If the variation among different time series is small, their 
predictive value is large. However, if the state changes of the individual time series 
are large and frequent, they are normally not representative of a system's behavior. 
Note that, by means of the time-dependent distribution of states (i.e. by simulation 
of the master equation alone), we cannot always distinguish between a strongly 
varying (e.g. oscillatory) system behavior and a "smooth" one. 

2. One can often distinguish desirable and undesired (maybe catastrophic) system 
states. In such cases we may be interested in the occurence probability with which 
the system under consideration up to time t takes a desired time series (containing 
desirable states only). This is relevant for prognoses as well as for the controllability 
of technical systems. 

3. Sometimes one also likes to know the occurence time distribution of a certain se- 
quencing of states. 

4. When occupying an undesired state, it may be interesting to evaluate the expected 
escape time until the system reaches, for the first time, one of the desired states. 

All these quantities can be obtained by evaluating a large number of stochastic time series. 
These can be generated and investigated in microsimulations applying the Monte-Carlo 
simulation method (cf. Binder 1979: Chap. 1.3.1). 

The Monte-Carlo technique utilizes the fact that a system under consideration stays 
in an occupied state z n _i for a time period r n — T n _i which is exponentially distributed 
according to (|3~4]). Therefore, in Monte-Carlo simulations the respective life time r n — T n _i 
is determined via the formula 

t„ - r n _i := — - — r In y n , (36) 

w(i n -i) 

where y n G [0, 1] are uniformly distributed random numbers (cf. Press et al. 1992: Chap. 
7). After this time period (i.e. at time r n ) the system goes over into one of the other states 
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i n 7^ i n -i with probability 

/. ,. \ w{i n \i n -i) P(i n ,r n + Ar|i n _i,r ft ) , . 

p(*n|»n-i) := — = Jim — • (37) 

W U\ i n-l) P(j,r n + AT\i n -i,T n ) 

3=1 3=1 
W/i»-l) (iV»n-l) 

The realized state i n is again selected by means of a uniformly distributed random variable 
z n G [0, 1]: The randomly resulting value z n corresponds to an occupation of the state i n 
for which 

In 

Z(z n _i) < z n < Z(i n ) with Z{i n ) = P{j<i n \i n -x):=^p{j\i n - 1 ) (38) 

3=1 

is fulfilled. The initial state io is usually given, but alternatively it can also be randomly 
chosen in a similar way. 

Time series which are generated by means of the above outlined method have the 
meaning of possible realizations of the stochastic process under consideration. They can 
be investigated in different ways. For illustrative reasons, let us assume to have generated 
A different time series Y a (r). If rij is the number of systems found in state j at time r, 
the distribution of states corresponds to 

n-(r) N 
P(j,r) = Jim^p since Y t n j {r)=A. (39) 

i=i 

For finite A, the probability P{n\, . . . , n^', t) to find the occupation numbers rij at time r 
is given by the multinomial distribution 

A\ N 

P(n a , ...,n N ;r) = : I] (40) 

ni! ...n N \ 

(cf. Helbing 1995: 70), so that the distribution of states nj{r)/A obtained by microsim- 
ulations will normally differ from P(j,r). This shows that we would need a tremendous 
number A of simulation runs to determine P(j, r) from microsimulations. For other quan- 
tities the situation is similar. As a consequence, the Monte-Carlo technique is rather a 
brute force than an efficient method. It is, however, very suitable for generating some sam- 
ple time series for illustrative purposes (which is sufficient for Item 1). The big advantage 
of the microsimulation technique is that it usually allows a relatively simple treatment 
of systems with a huge number of possible system states and/or very complex relations 
for the transition rates. Therefore, it is mainly used in situations where the derivation or 
simulation of a master equation is too difficult. 

Since the microsimulation technique normally requires an extreme computational ef- 
fort, it should only be applied if there are no suitable alternatives. As discussed in Section 



372|, the distribution of states P(j, r) is much easier obtained by simulating the master 
equation. In the next section it will be shown that more efficient methods than the mi- 
crosimulation technique can also be developed for the problems raised in Items 2, 3, and 
4. The reason is, that they only regard the sequencing of time series up to a certain time 
t but not the waiting times after which the single state of the time series are occupied. 
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4.2 Occurence Probabilities and Occurence Times of Paths 

With the results of sections [| and |3] we can now calculate the probability P(H n ,t) that 
we have the event history 

H n := {t , io; h, ii; . . . ; t n , i n } = {t"o, «o; n, i\, . . . ; r n , z n } (41) 

at time r(i). Remembering that S k (Tk\ik-i) — P{Tk > r k\^k-\ — h-i) with 

% ■= r(T k ) (42) 

is the probability to stay in state up to time r k and that ii;(zfc|zfc_i)<iTfc = P(ik,Tk + 
d,Tk\ik-i, Tfc) is the probability to change from state i^-i to state Z& between times and 
r k + <ir fc , we find 

P(H n ,r) = 5 n+1 (r|i ri )u;(2 n |i n _i)(ir n S'' 1 (r n |z ri _i)u;(ii n „i|z„_2)^n-i 

. . .w(i 2 \i 1 )dT 2 S 2 (T 2 \ii)w(ii\i )dT 1 S 1 (Ti\i )P{i , To) (43) 

(Empacher 1992). Often, however, one is not interested in the times = r(i&) at which 
the single transitions occur, but only in the sequencing, i.e. in the path 

C n ■= in < in-1 < • ■ ■ < H< ^0 (44) 

which the system has taken up to time r(t). For example, one could ask which is the 
probability that somebody is married the second time or unemployed the third time at 
time r. One could also compare the probabilities of being married the second time after 
divorce or after death of the partner. The necessary quantities for answering questions 
like these can be derived from formula fl4"3f) by integration with respect to r 1; . . . , r n . For 
the probability of having path C n at time r we obtain 

P(C n ,r) = J2 — n Pm k (wk,r)w(C n )P(io,t ) . (45) 

fc =° J] {wi-w k ) 

l=o 

Here, is the multiplicity of the overall departure rate Wk : = w(ik) in path C n , 

f 1 if n = 

w(C n ) = w(i n <- i n -! <— . . . <- z ) := 



(46) 



11^(^1*1-1) if n > 1, 

1=1 

/_i\m+l / m-2 (m-l-m) / nj-1 „(m-n 2 ) 

^ = K^+E- K w) +E- 

m(m-l)\ ^ ni \ nt^i n 2 



for m > 2, p m = 1 for m = 1, and 



/ +1 W):= E 7 vTT-^o (4* 
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with 5io — 1 if I — and 810 = otherwise. The detailed steps on the way of deriving 
P(C n ,r) are presented in a paper by Helbing and Molini (1995). If one restricts to the 
case that all overall departure rates Wk are different from each other (i.e. to pure birth 
processes) the above relation simplifies to 

P(C n , t) = £ — w(C n )P(t , r ) , where S 1 ^) = e^^o) . 

fc=0 II Mil) - w (ik)\ 

1=0 

CM*) 

(49) 

The special formula (HS) was already presented by Chiang (1968: 50ff) who, however, did 
not introduce the more general path concept developed above. 

The probability density P(r\C n ) of the occurence times r (the occurence time distribu- 
tion) of path C n can now easily be obtained from (|45|). It is given by the formula 

P(C t) 7 n 1 

P(r\C n )= with P(C n ):= dTP(C n ,T)=w(C n )P(i ,T )H—-- (50) 

(Helbing, 1994, 1995). From this we can derive the average 

(T ) c „ = / dTT P( T | C „) = £ (51 ) 



TO 



fc=0 



o/ t/ie occurence times r and their variance 



n i 

(cf. Helbing 1994, 1995). The average occurence times are the basis for the calculation of 
the expected escape times. 



4.3 The Simulation Tool EPIS 



For a determination of the quantities called for by Items 2 and 4 of Section [0] we have 
to sum up over the occurence probabilities or average occurence times of many paths. In 
the following, we will discuss how this can be done in an efficient way with respect to 
computer memory and simulation time. For illustrative reasons we begin with the socalled 
contracted path integral solution of the master equation. This reads 



oo 

^» = EE*r) (53) 

n=0 C n 



with 

N N N 

£:= E E ••• £ > (54) 

C n 8 n-l = 1 8 n-2 = 1 l = 1 

(i n -l&) (»n-27^n-l) ('05**1) 

since, using r = 0, the probability P(J, r) to find state j at time r is given as the sum 
over the occurence probabilities P(C n , r) of all paths C n with an arbitrary length n which 
lead to state i n := j within the time interval r (Helbing 1994, 1995). 
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When ( |53|) is evaluated numerically, one must restrict the summation to a finite number 
of relevant paths. Here, we can utilize the fact that the occurence probability P(C n , r) of 
a path C n with a non-absorbing final state j is negligible if 

\t - (r) Cn \ < ayfoZ (55) 

with a suitably chosen accuracy parameter a. For r < (r)c n — cl\J ©c„ there is not enough 

time to traverse all states ik of the path, whereas for r > (r)c n + a J ©c n the system 
will probably have already visited additional states i n +i, i n +2, ■ ■ ■ Therefore, only paths 
which fulfil condition ( [55]) are relevant. Choosing the accuracy parameter a « 3 allows 
to reconstruct about 99% of the probability distribution P(j, r) which can be checked by 
means of the normalization condition Y,f=i P{jj t) = I. 

If we are interested in the probability that the system takes desired paths only, we 
have to take into account an even smaller number of relevant paths. This is simply done 
by restricting the summation (|4]) to the desired states. A similar procedure is used to 
determine the expected escape time of a system from undesired states. 

In order to evaluate the path-related quantities introduced above, the simulation tool 
EPIS (Efficient Path Integral Simulator) has recently been developed at the University 
of Stuttgart (Molini 1995). This bases on the standard path-search algorithm depth-first 
(Schildt 1990) but restricts it to the relevant paths (see Figure 2). The advantage of this 
algorithm is that it uses computer memory very efficiently and allows to calculate certain 
quantities like 

( r )c n+1 = (r)cn + — , @c„+i = @c„ + ( 1 s 2 , w{C n+l ) = w{C n )w{i n+1 \i n ) (56) 

after each step. In this way, the results of previous calculations can be utilized which again 
saves computer time. 

The simulation tool EPIS is constructed for systematically generating relevant paths 
and for calculating 

• their occurence probabilities (fOp, 

• their occurence time distributions (|50|), 

• the occurence probability of desired paths up to time r, 

• the expected escape time from undesired states. 

Moreover, the contracted path integral solution ( |53| ) can be compared with the distribution 
P(j, t) of states which results by a numerical integration of the master equation. 

5 Effective Cumulative Life-Time Distribution and 
Hidden State Concept 

The time-dependence of hazard rates is usually fitted to functions like the ones be- 
longing to the Gompertz distribution, the Weibull distribution, the log-normal distribu- 
tion, the log-logistic distribution, the sickle distribution, the extreme-value distribution, 
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Step 1 



Steps 2-4 



Steps 5-11 



Steps 12-14 



Figure 2: Efficient path generation according to the modified depth-first algorithm. Here, it is 
assumed that the system can be in one of the states Y n E M = {1, 2, 3} after the nth transition 
and that it is initially in the state Yq = iq = 2. The single steps of the procedure are symbolized 
by arrows. They correspond to an extension of a path C n by a new admissible state Y n+ \ = i n+ i 
or to a removal of the last state i n if the path C n cannot be further extended to a relevant 
path C n+ \. Full circles represent the states of the last generated path and are connected by 
thick arrows or lines. Crossed circles (®) stand for states Y n+ \ which are not admissible since 
transitions from state Y n must lead to another state Y n+ \ ^ Y n . The symbol "0" indicates that 
the resulting path C n +i turns out to be not relevant due to r > (r)c n+1 + ay / 0c n+1 . Since then 
all longer paths are irrelevant, too, the procedure removes the last state and tries to extend the 
remaining path C n by another state. If this is not possible, the state i n is also removed, etc. 

the gamma distribution, and others (Blossfeld et al. 1986; Kalbfleisch/Prentice 1980; 
Elandt- Johnson/ Johnson 1980; Tuma/Hannan 1984; Cox/Oakes 1984; Lancaster 1990; 
Courgeau/Lelievre 1992). Since these functional time-dependences can not always be in- 
terpreted in terms of underlying social or psychological processes, a new concept for an 
interpretation of the respective time-dependence of transition rates is proposed in the 
following. This bases on the discovery that, if we would have sequential or parallel tran- 
sitions between hidden states, this would cause a time- dependence of the related effective 
transition rate w e s(r;j\i) of the total transition process. Hidden states denote states which 
are either not phenomenologically distinguishable from other states, not externally mea- 
surable, or simply not detected. They could, for example, reflect the psychological and 
mental stages preceding a concrete action, the individual predisposition, or the internal 
attitude towards a certain behavior. 

It should be underlined that the aim of the hidden state concept is not to question 
or criticize but to supplement the powerful methods of classical survival analysis, which 
have proved their suitability in numerous applications. Favourable properties of the hidden 
state concept are 

1. the huge number of classes of effective transition rates which can be generated by 
it and which may include additional kinds of time-dependences, 

2. the theoretical relationship between these different classes of time-dependent tran- 
sition rates, 

3. the possibility of a direct interpretation of each hidden state model. 
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Before the hidden state concept is discussed in general (cf. Section |573|) , it will be 
illustrated by two well-known cases which correspond to pure birth processes (cf. Section 



PD and to the situation of unobserved heterogeneity (cf. Section |5.2| ), respectively. This 
allows to introduce the basic ideas without complex mathematical considerations which 
are unavoidable in the general case. 

5.1 Sequential Transitions 

As an example for a series of birth processes, let us face the transition 

married — > divorced (57) 

and assume that it is actually a transition of the kind 

married — ► married with attitude towards being divorced — ► divorced . (58) 

It should be mentioned that this example is, of course, only a "toy model" which was 
chosen exclusively for illustrative reasons. 

In order to take into account the hidden state "married with attitude towards being 
divorced" we can introduce the following notation of states: 



(1,1) 
(1,2) 
(2,1) 



married , 

married with attitude towards being divorced , 

divorced . (59) 



With this, example ( |58|) can be written in the form 

H 2 :=(2,1)<— (1,2) <— (1,1) (60) 



or 



T-(-2-=}< — ii < — io, (61) 

where j = h ■= (2,1), ii := (1,2), and io := (1,1). In general, the distinguished states 
% G {1, . . . , iV} actually comprise the "hidden states" 

i:=(i,h) with h G {1, . . . , N*} , (62) 

where h differentiates the (possibly unknown) variants of state i. 

Analogous to ( |P2"D the effective cumulative life-time distribution F^(r\i) for the total 
sequence fl6l]) of transitions can be defined by the probability P{T 2 < t\\\ <— io) that the 
final transition (i.e. the transition into state j) takes place before time r given that the 
preceding path was ii « — i . Since P(TC 2 ,r) is the probability to have path H 2 at time 
r we get 

FMi) = P(T 2 < r\h <- i ) = P{H 2 , r) . (63) 



Inserting ( f49|) leads to 

( S\r\i Q ) S\r\i x ) 



-w(i )[w(ii) - iw(io)] -tw(i 1 )[u;(io) - «;(ii)] 

+ T. ! T . A w(h)w(i ) (64) 
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because of u>(j) = 0, = w(ii), w(ii|io) = w(h), an d P(h, T o) — 1- With this we 

can define the effective survivor function 



SMi) := 1 - Fi(r|0 = f ^ - ^) : (65) 

V ^( J o) / w(ii) -w(io) 



in accordance with 

In the following we try to represent this expression in the form 



offl r l*J = CX P 



(66) 



(cf. flTDp) which corresponds to the detected (observed, measured) transition ( E7| ) of the 
type J < — i (where i corresponds to "married" and j to "divorced"). The appropriate 
relation for the effective overall departure rate w e s(r\i) is found by differentiation of (|66|) 
with respect to r. We obtain 

-^Sls(r\i) = -w e s(r\i)Sl e (T\i) i.e. w eS (r\i) = -^y^y ^^(^1*) • ( 67 ) 

Inserting ( |65"1) finally gives 

w os(^) = Bunri — "eu'Si'\ ■ ( 68 ) 



s 1 


(r|io) - 


S 1 


V|ii) 


s 1 


>|io) 


S 1 





w(i ) w(ii) 

For the effective transition rate w e s(r; we find w c s(r;j\i) = w e s(r\i) since there are 
no alternative transitions from state % to state j. 

The most important fact about the relations for w e ff(r;j\i) and w e g(r\i) is that these 
are time- dependent. Therefore, the hidden state concept of behavioral changes may serve 
as a means for interpreting the time-dependence of hazard rates. A comparison with em- 
pirical data, however, shows that the above model fl58|) is still too simple for an explanation 
of the sickle curve which is found for the hazard rate of divorces (Diekmann/Mitter 1984). 

The above toy model could be improved by distinguishing different motivations (i.e. 
alternative reasons) which may lead to a divorce. This brings us to hidden state models 
with parallel transitions which correspond to cases of unobserved heterogeneity. 



5.2 Parallel Transitions 

For illustrative reasons we will discuss an example which stems from Kinsey et al. (1948) 
concerning the kinds of sexual activities of white males in the United States. It should be 
mentioned that this rather old example is not intended to demonstrate survival analytical 
or statistical methods. In addition, since it is based on a small amount of cross-sectional 
data, it only allows the determination of the most frequent transitions. Longitudinal 
data would, of course, supply more detailed information, e.g. concerning the possible 
(but seemingly infrequent) transitions "heterosexual activity" — ► "bisexual activity" or 
"bisexual activity" — > "homosexual activity". Therefore, the interpretation of the data 
will also not be the focus of our discussion. The example was rather chosen for didactical 
reasons, since it is suitable for illuminating 
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1. the hidden state model with exclusively parallel transitions, 

2. possible methods for the detection and separation of hidden variables (here: two 
different kinds of "homosexual activity" ) from empirical data, 

3. a simulation model basing on the master equation. 



100 




20 25 
Age (years) 

Figure 3: Percent of white males engaging in heterosexual activity (□), bisexual activity (A), 
homosexual activity (O), or having no socio-sexual contacts (x) in dependence of age. Here, 
"bisexual activity" summarizes the five classes "predominantly heterosexual, only incidentally 
homosexual" , "predominantly heterosexual, more than incidentally homosexual" , "equally het- 
erosexual and homosexual" , "predominantly homosexual, more than incidentally heterosexual" , 
and "predominantly homosexual, only incidentally heterosexual" . The solid lines show the sim- 
ulation results of the hidden state model which is proposed later on. Despite its simplicity it 
apparently fits the data quite well. 

We now come to the example by Kinsey et al. (1948). Figure 3 shows the distribution 
of the different kinds of sexual activities (within a period of three years) in dependence 
of age r. It turns out that the proportion 



X tot (r) :=£X,(r) 
j'=i 



(69) 



of males with socio-sexual contacts can be approximated by the logistic equation (Pearl 
1924; Verhulst 1845) 

^ = uX tot (r)[l - X tot (r)] , (70) 
where X tot (0 years) = 0.0173, v = 0.36 per year, and 

Xi(t) := proportion of white males engaging in heterosexual activity (j — 1) , 
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X%(t) : = proportion of white males engaging in bisexual activity (j = 2) , 

X 3 (t) := proportion of white males engaging in homosexual activity (j = 3) . (71) 



Equation (|70|) could be interpreted in the way that the proportion 1— X to t of males without 
socio-sexual contacts is seduced to sexual contacts by the proportion X tot of individuals 
of about the same age with socio-sexual contacts. 

However, note that the logistic curve deviates from the empirical data for very young 
and for old males. Whereas the proportion of males without socio-sexual contacts should 
start with a value of 1 at the age of years, it should increase after an age of about 30 
years due to the decrease of sexual opportunities (or other reasons). In any case, empirical 
and theoretical statements are very questionable before an age of 10 or 15 years. 

We will now investigate the proportions 

P (j,r):=^\ (72) 

of white males with socio-sexual contacts who engage in sexual activity j. These are 
represented half-logarithmically in Figure 4. P(j, r) can also be interpreted as probability 
that, in a representative sample of males with socio-sexual contacts, a randomly picked 
out male engages in sexual activity j. 




l 1 1 1 ' 1 1 1 1 

10 15 20 25 30 35 40 45 

Age (years) 

Figure 4: Half-logarithmic plot of the proportions of sexually active males engaging in hetero- 
sexual activity (□), bisexual activity (A), and homosexual activity (O). The broken lines show 
the corresponding linear regression curves (on the basis of the time period 10 years < r < 45 
years for bisexual activity and the time period 25 years < r < 45 years otherwise). The numbers 
in brackets are the respective correlation coefficients. 

We find that heterosexual activity increases, but bisexual and homosexual activity 
decrease with advancing age which might be a consequence of an adaptation to social 
norms. The half-logarithmically scaled curve of bisexual activity is quite well decribed by 
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a linear relation over the whole range of ages, whereas the other curves behave almost 
linearly after an age of 25 years. Between 10 and 25 years there seems to be a surplus of 
homosexual activity connected with less heterosexual activity than expected. This surplus 
decays very fast and fills the gap of heterosexual activity. We may therefore suspect to be 
confronted with two different types of homosexual activity: 

1. a rather permanent type of homosexual activity (with a half-life period of about 30 
years) which might arise from a "homosexual predisposition" , 

2. a short-lived type of "adolescent homosexual activity" (with a half-life period of 
about 2 years) which might be a substitute for lacking opportunities for heterosexual 
activity. 

This interpretation is in good agreement with the findings of developmental psychology. 

Summarizing the previous considerations we predominantly have the following transi- 
tions 




between the states 



j 


:= (1,1 




:=(2,1 




:=(3,1 


3 


:=(3,2 



= "heterosexual activity" , 

= "bisexual activity" , (73) 

= "homosexual activity with homosexual predisposition" , 

= "adolescent homosexual activity, but heterosexual predisposition" , 



where we have again used the hidden state representation flB2"|). As Figure 3 shows, the 
temporal course of the proportions Xj(r) can, for r > r := 10 years, be approximated 
by ( [70"D and ( [72] ) together with the master equation 

^p(i', t) = £ [Mi'|i)P(i, t) - Mi|i')P(i', t)] (i, i' g {j, il, il, ijj}) . (74) 

The presented simulation results are for the parameter values 

P(iJ,T ) = 0.33, = 0.042 per year, 

P(gr ) = 0.05, = 0.024 per year, 

P(ig, r ) = 0.35 , = 0.366 per year (75) 

which were determined with the method of least squares (Elandt- Johnson/ Johnson 1980; 
Tuma/Hannan 1984; Lancaster 1990; Helbing 1995). Although the terms to(i'|i)P(i, t) 
were assumed to be negligible for i' ^ j in accordance with the simplified model under 
consideration, the simulation results fit the empirical data quite well. 

A number of alternative models with the same number of parameters have also been 
tested. Most of them yielded worse correlations with the empirical data. Only the assump- 
tion that "homosexual activity with homosexual predisposition" is replaced by "bisexual 
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activity" instead of "heterosexual activity" produced comparable results. A better cor- 
relation with the empirical data can, of course, always be reached by more complicated 
models which include a greater number of parameters. However, an increase of the number 
of parameters is only advisable, if this gives a considerably better fit. Otherwise, some 
of the parameters will be insignificant so that the explanatory power of at least some 
variables is very low. 

In the following we want to describe the transitions from the two "hidden" homosex- 
ual states 1q = (3, 1) and ijj = (3, 2) in an overall manner as if we would only have one 
homosexual state i — 3. The associated effective transition scheme (which corresponds to 
the states that were actually detected) is 

2 

:x i 

3 

The effective cumulative life-time distribution F^ s (r\i) = P(Ti < t\X = i) for the transi- 
tion from homosexual to heterosexual activity is given by the probability that one of the 
paths j < — ip or j < — ijj is taken before time r: 

Fi(r|<) = P(j - i* r) + P(j «- ig, r) . (76) 

From ( pE9"| ) we obtain 

= [l-S\r\il)]P(gT ) + [l-S\r\i 3 )]P^,r ) (77) 

due to = w(il), w(j|ig) = u;(ig), and w(j) = 0. With 5^ ff (r|z) = 1 - F^(r\i) we 

now find the effective survivor function 

SlMi) = S\r\il)P(il, r ) + S\r\il)P(il, r ) . (78) 

If we again try to express the effective survivor function in form (|B"6| ), we can calculate 
the effective overall departure rate w e s(r\i) via formula (|6~7|). Inserting (ff§|) provides 



w(il)S\r\ 


il)P(f ,r )+w(il)S\r\ 






i§)P(i§,7b)+5*(r| 


'4)P( 


'M To) 



This result obviously differs from expression fl68"|) which we obtained for the hidden state 
model ( pT| ) of the Section |5.1| . However, since we have only transitions to one final state, 
we again find w e s(r;j\i) = w e s(r\i) for the effective transition rates w e s(r; 

Note that the two hidden homosexual states are an example for the widespread case 
of unobserved heterogeneity (unmeasured heterogeneity) (Blossfeld et al. 1986; Elandt- 
Johnson/ Johnson 1980; Diekmann/Mitter 1984; Tuma/Hannan 1984; Cox/Oakes 1984; 
Lancaster 1990; Courgeau/Lelievre 1992). In the general case of heterogeneous populations 
(mixed populations), we have the relations 

S 1 cS (r^)=T,S 1 (r\i k )P(i k ,r ) (80) 

k 
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and 

E/^JIioWoVo) 
Wesir; j\i) = w e s(T\i) = — ; ; (81) 

CSV ,J\ ) efn I , <£ S \r\i*)P(%,T ) ^ ' 

k 

if there is only one final state j = (J, 1) (Courgeau/Lelievre 1992: 44; Blossfeld et al. 1986: 
93). Formula ([80]) describes the mixture of survival functions in heterogeneous popula- 
tions. P(iq,Tq) is called the mixing distribution (compounding distribution), and the state 
i — iq summarizing the states ig can be interpreted as compound state. In the special 
case that some unobserved subpopulations (which are here reflected by k) are subject 
to transitions but others not, we are confronted with a so-called mover-stayer model (cf. 
Courgeau/Lelievre 1992). The reason is that the formerly mentioned subpopulations can 
be interpreted as movers, the latter ones (which are characterized by w(iq) = 0) as stayers. 



5.3 The General Case 

An investigation of formula d68|) shows that sequential transitions are related with an 
effective transition rate w e g(r;j\i) which is zero at r = r and increases monotonically in 
the course of time up to the limiting value min[w(i ), w(ii)]. In contrast, for the case (|79|) 
of parallel transitions we find that the corresponding effective transition rate w e g(T;j\i) 
starts with the finite value w e g(r ]j\i) = w(j^)P(t^, r ) + u;(ig)P(ig, r ) and decreases 
monotonically up to the lower limit min[w(io), w (ig)]. 

Therefore one can conjecture that any time-dependence of hazard rates (including 
sickle-shaped ones) can be obtained or at least approximated by a suitable combination 
of sequential and parallel transitions. The general hidden state concept being necessary 
for this is presented in the following and can be skipped by readers who are not interested 
in the mathematical details. 

Again we assume that the detected states i comprise a number of hidden states i = 
(i, h) with h G {1,2,..., Ni} and Ni > 1. In the general case, a detected transition j < — i 
is composed of all kinds of transitions between the hidden states ifc = (i, hk) comprised 
by state i (hk G {1, . . . , Ni}) before the individual or system comes the first time into one 
of the hidden states j = (j, h n ) (h n G {1, . . . , Nj}) which are comprised by state j (see 
Figure 5). Therefore, we will have to sum up over all paths TC of the form 

H n :=j< — i n _i < — . . . < — i with i k = (i, h k ) and j = (j, h n ) , (82) 

where the path length n is arbitrary. For this purpose we introduce the abbreviation 

oo Nj Ni Ni Ni 

j:--=j:j:j: ^ e== e e ••• e (83) 

H n=0h n = lH n H n h n - 1 = l h»-2=l h =1 

(fc„_2=?^ n _l) Chit^i) 

similar to (|54]). 

We are now looking for an expression which allows to calculate the probability 
F^(t; j\i) that the system changes, for the first time, to the detected state j before 
time r given that the preceding detected state was i ^ j (cf. (|i~3D). Here, we remember 
that P(TC n , t) is the probability that the path 7i n is traversed before time r, but not 
extended by additional states (cf. (f49D). Since setting w(j) = guarantees that the final 
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Figure 5: General hidden state model including parallel and sequential transitions illustrated for 
the case iVj = 3 and Nj = 2. The possible hidden transitions leading from the detected state 
i £ {(i, 1), (i, 2), (i, 3)} to the detected state j € {(j, 1), (j, 2)} are represented by arrows. The 
thick arrows indicate one of the numerous paths Ji n of the form ( |82|) which the system may 
take. 



state j of H n is not left any more, P(H n ,T\w(j) = 0) is the probability that the path 
7i n is traversed before time r. For this reason, F^ s (T;j\i) is given by the sum over the 
probabilities P(TC n , T \w(j) — 0) of all paths 7i n which have the form 



oo Nj 



F&t; = £ P(n n , r\w(i) = 0) = E E E r\w(i) = 0) . (84) 

H r n=0 h n =l H n 

The factor P(i) takes into account that F^ s (r;j\i) is a conditional probability. It is deter- 
mined by the probability 

N oo Nj 

P(i) :=EEE Eig,*^(i) = °) ( 85 ) 

3=1 n=0 h„=l H n 

that the system starts with the detected state % and changes to another state j ^ i at all. 
It is obvious that the evaluation of formula (|8~4]) requires a computer program like EPIS. 
From relation (|54"D we can obtain the effective cumulative life-time distribution 

N 

F«{t\{):= J2F^(r;j\i) (86) 

3=1 

analogous to ( |I~4"D and confirm the desired normalization condition 

JkaFi(T|i) = l. (87) 
Moreover, we can derive the effective survivor function 

Sl s (r\i) := 1 - Fi(r|i) , (88) 
the effective •probability density function 

fe S (r;j\i)-=^ s (r;j\i), (89) 
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and the effective transition rates 

w eS (r;j\i) == / S ( ?n ) (90) 

in accordance with ( [16D and (|T7|). For the effective overall departure rates 

N 

w e s{r\i) := E w eS (r;j\i) (91) 



i=i 



we again find the simple relation 

In their general form, the above formulas cannot be further evaluated. For illustrative 
reasons we will calculate the effective overall departure rates for the special case that there 
is exactly one hidden path of the form fl52|) which corresponds to the detected transition 
j < — i. This implies that the length n of the path as well as the states and j are 
uniquely determined by i and j, i.e. n = n(i,j) and h k = h k (i,j) (k G {1,2, ... ,n}). 
Therefore, the hidden states ii,...,i n _i could be interpreted as transient states. If all 
overall departure rates w k = w(i k ) are pairwise different from each other, we find 

P(H n ,T\w(j) = 0) 

/ 71-1 ^(tW, ) 1 \ 

= E + ] ^(Wn)P(io, ro) (92) 

v k =° [o - w(i k )] n [w(h) - w(i k )] n - o] 

and, due to w(i k ) > 0, 



1=0 l =o 




lim P(H n , r\w{j) = 0) = — ^ w{H n )P{i , r ) . (93) 

n 

i=0 



+ — E w(H n )P(i ,r ) 



n-l 

fc= ° -^(ifc) n H^) ~ ^( i fe)i n i °( i o ' 

F&r\i) = ^ a (94) 

— w(n n )P(i ,r ) 

Yl w(i{) yyi) 

1=0 

(which can be further simplified) and 

S\r\i k ) 



/n-l \ n-l 

Sl s (r\i) = wfr) E 

\i=0 / fc=0 



n-l 

w(i fc ) - itf(ifc) 



(95) 



i=0 

(i^fc) 
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Slg(r\i) is obviously a direct generalization of S 1 (r\i) since, for n — 1, we find 

5i(r|i) = ^(io)^^ = ^(r|i ) = S\r\i) . (96) 
Finally, we arrive at the desired result 

n-1 



fc=0 



n m^) - 



Weff(rli) = — ; (97) 

^(tIU) 1 ^ 



E 

fc=0 



n-1 



**) II [w(k) - w(i k )} 



1=0 



due to 

^-S 1 (r\i k ) = -w(i k )S 1 (r\i k ). (98) 

Although the special formula ( j97[) is restricted to hidden state models consisting of se- 
quences of n birth processes, it is already much more complicated than ( |68"D . 



6 Summary and Outlook 

In this paper we showed that, in the Markov case, survival analysis is related with the 
master equation which describes the temporal evolution of the distribution of states. The 
simulation of the master equation is a powerful technique for the investigation of various 
stochastically behaving systems in physics, chemistry, biology, and the social sciences. It 
is particularly suited for scenario techniques, since the numerical integration of the master 
equation is normally quite simple and fast. 

However, some interesting questions related with the longitudinal time series of 
stochastic processes cannot be answered by means of the master equation. Whereas the 
microsimulation of sample time series can be done with the Monte-Carlo technique, it is 
rather inefficient for the numerical determination of quantities related with the sequenc- 
ing of time series. Therefore, the simulation tool EPIS has recently been developed at the 
University of Stuttgart. It facilitates the generation of relevant paths and the evaluation 
of formulas which we were able to derive for path-related quantities. This includes the 
occurence probabilities and occurence time distributions of paths, or the expected escape 
time from undesired states. 

Finally, the formula for the occurence probabilities of paths allowed to develop a 
hidden state concept of behavioral changes which can serve as a means for interpreting the 
respective time-dependence of hazard rates. Starting from a certain hidden state model, 
it is possible to derive the corresponding effective transition rates which can be compared 
with the time-dependence of the empirically obtained hazard rates. The different steps 
which are necessary for determining a suitable hidden state model and the corresponding 
parameter values (including the hazard rates) were illustrated by a concrete example. 

Present research focuses on the investigation of the following questions: 
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1. Which kinds of time-dependences can be interpreted or approximated in terms of a 
hidden state model? 

2. Does the time-dependence of hazard rates determine the corresponding hidden state 
model in a unique way? 

3. If not, which are the transition schemes of the alternative hidden state models and 
which of them is the simplest or most plausible one in terms of a sociological or 
psychological interpretation? 
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